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ABSTRACT 

A two dimensional hydrochemical hybrid code, KM2, is constructed to deal with astrophysical 
problems that would require coupled hydrodynamical and chemical evolution. The code assumes 
axisymmetry in cylindrical coordinate system, and consists of two modules: a hydrodynamics module 
and a chemistry module. The hydrodynamics module solves hydrodynamics using a Godunov-type 
finite volume scheme and treats included chemical species as passively advected scalars. The chemistry 
module implicitly solves non-equilibrium chemistry and change of the energy due to thermal processes 
with transfer of external ultraviolet radiation. Self-shielding effects on photodissociation of CO and 
H 2 are included. In this introductory paper, the adopted numerical method is presented, along 
with code verifications using the hydrodynamics module, and a benchmark on the chemistry module 
with reactions specific to a photon-dominated region (PDR). Finally, as an example of the expected 
capability, the hydrochemical evolution of a PDR is presented based on the PDR benchmark. 

Subject headings: Astrochemistry — Hydrodynamics — Methods: numerical — Photon-dominated 
region 


1. INTRODUCTION 

Spectral lines from molecules or ions in the interstellar 
medium can help to probe physical properties of the gas 
through density, temperature and kinematics. Chemical 
modeling of the interstellar medium is crucial to have 
accurate interpretations of astronomical observations as 
the intensities of spectral lines depend on spatial distri¬ 
bution and abundances of the individual species. Recent 
advances in computing power also enable comprehensive 
radiative transfer simulations to predict a large number 
of observable s pectral lines from theoretical and chemi¬ 
cal rnodels ('e.g..lJuvelalI99^lHogerheiide fc van der TakI 
l200r)HBrinch fc Hogerheiideil20ir)ll . Next generation tele- 

scopes such as the Atacama Millimeter Array (ALMA) 
will significantly improve the sensitivity, angular and 
spectral resolutions of the observations. It is essential 
to establish accurate physical and chemical models of 
the interstellar media to derive physical properties of ob¬ 
served objects by comparing observations and theoretical 
models. 

Chemical evolution of the interstellar medium is tightly 
coupled with its hydrodynamical evolution. The rates of 
chemical reactions depend sensitively on the density and 
temperature of the gas. Molecular cloud cores with dif¬ 
ferent magnetic field strengths, angular momenta, and 
initial density profiles have different hydrodynamical 
evolutions during gravitational collapse. Consequently, 
these cores are expected to have different chemical evo¬ 
lution as a result of the different distributions of den¬ 
sity and temperature. Recent observations reveal a va¬ 
riety of chemical properties in star forming regions that 
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might have hinted such connections. For instance, de¬ 
pletion of CO is found at the dense interior of many 
starless cores, while L1521E does not show sign of de¬ 
pletion in spite of its central de nsity similar to other 
objects (jTafalla fc Santiagoll2004fl . Some protostar s har¬ 
bor complex organic molecules (iCazaux et al.ll200^. oth - 
ers harbor carbon-chain molecules ( Sakai et al.l I200§1 . 


These observed variations of chemical properties might 
be caused by their respective hydrodynamical evolution 
inside each of the cores. Hydrochemical modeling of star 
formation process may help with the insights into the 
origins of the observed differences. 

Many simulation works have been done to investigate 
chem i cal evolution in star forming cores (lAikawa et aH 
l2005t Ivan Weeren et al.l 120091 Fnruva et al.ll20I2tl . cir- 
cumstellar disks ( Visser et al.l 120091 1201 III , and so on. 
In these studies, chemical evolution was solved indepen¬ 
dently of hydrodynamical evolution using the density, 
temperature, and velocity obtained by hydrodynamic 
simulation or semi-analytic model. Therefore, chemi¬ 
cal evolution did not affect thermal and hydrodynami¬ 
cal evolutions in their simulations. However, chemistry 
plays an important role in some astrophysical flows, be¬ 
cause chemical abundance of the gas affects thermody¬ 
namics, and consequently dynamics, of flows. For in¬ 
stance, atomic or molecular cooling is cr ucial for shock 
induc ed formati on of molecular clouds (lAsahina et al.l 
I20I4I1 and stars (jVanhala fc CameronllI99^ . As shown 
in this paper, photoevaporation flow of interstellar gas 
caused by strong external radiation field is also typical 
example in which chemical evolution affects thermal and 
hydrodynamical evolutions. In these situations, chem¬ 
istry and hydrodynamics should be solved together in 
self-consistent way. 

Most previous works with hydrodynamic simulations 
including chemical reactions have focused on evolution 
of primordial gas in the cosmological context. Chem¬ 
istry is relatively simple in the primordial gas as it only 
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involves about ten species consisting of primordial ele¬ 
ments such as hydrogen, helium, and deuterium. On 
the contrary, chemical networks in the present-day inter¬ 
stellar medium are larger and more complex. To date, 
more than 170 interstellar molecules have been listed 
in the Cologne Database for Molecular Spectroscopy 
(CDM^). Substantial efforts have been made to develop 
systematic solvers to follow the chemical evolutions of the 
species that would have played roles in the astrophysi- 
cal problems of intere sts. For example. Na hoon (part of 
the ” KID A” package (lWakela,m et al.ll2 01^. Wa kelam & 
Herbst 2008), and KROME ( Grass! et al.ll201^ are two 
such packages that have incorporated chemistry solvers 
to solve large networks of chemical reactions. Hydro- 
dynamic simulations aimed at the present-day universe 
would have to employ systematic solvers of this kind to 
properly tackle the problems. 

Hydrochemical simulations of present-day universe 
have to solve multi-physics characterized by various dif¬ 
ferent time scales, namely hydrodynamics, chemical re¬ 
actions, thermal processes, plus the radiative transfer. 
Hydrodynamics certainly affects chemical reactions in in¬ 
terstellar medium, while chemical reactions may play a 
role in the composition and thermal properties of the 
gas. Radiation may drive evolution of the interstellar 
medium through photoreactions and heating processes. 
Although integrating all the physics directly with respect 
to time by their smallest time scale is straightforward and 
a simple way to achieve numerical stability and reliable 
results, long-term evolution is prohibited. It is necessary 
to establish an effective time integration method main¬ 
tain both computational efficiency and accuracy. For 
these reasons, we developed a robust and computation¬ 
ally effective hydrochemical hybrid code, KM2, which we 
introduce in this paper. 

We will adopt models for the Photon-dominated 
regions (PDRs) as our benchmarks for capabili¬ 
ties developed with our hydrochemical hybrid KM2 
code for several reasons. First and foremost, the 
physical and chemical properties are dominated 
by penetrating far ultraviolet radiation (FUV: 
Q < hv < 13.6 eV), and their characteristics and 
emissions have been studied by a few code specif¬ 
ically developed f or th e m: Cloudy (iFerland et al] 

1199^ iShaw et all 120051: lAbel et al. 1200511 . Meudon 


{ Le^^ourlo^^^^ 1 19931 LePetit_e£jh 1200611. U CL_PDR 
( Papadopoulos et al.l I2002t IBell et al.l 20 05ll. Leiden 


I Black fc van Dishoectd TlOsW Ivan Dishoeck fc Black 

1988t iJansen et ahT 199511. CO S TAR ( Kamp fc Bertoldi 


Kamp fc van ZadelhoftI 1200111 . and 3D-PDR 
( Bisbas et al l 1201^ and so on. Comparisons and 
benchmark s tudies of these PDR codes were made by 
iRbllig et al.l (j2007tl using the PDR models based on 
a common list of reaction networks. Although the 
hydrodynamics was not solved in these codes tested 
for the benchmarks, it is useful to verify our solvers 
for chemical reactions, effect of radiative transfer, and 
thermal processes included. We first verified our code 
against these PDR benchmark models by performing 
chemical model calculations following similar criteria 
without actually advancing the hydrodynamics. Once 
the behaviors of the chemistry-only modules have 


http://www.astro.uni-koeln.de/cdms/molecules 


been confirmed, the hydrodynamics module have also 
been tested following a few classical examples. A 
photoevaporating PDR model in ID, with both of 
the hydrodynamics and chemistry modules active, was 
performed to showcase the potential dynamical effects 
from the hydrodynamics. To our knowledge, this is the 
first study of an evolving PDR model along with the 
impacts of its background thermal flow. 

This paper aims to present an efficient approach to 
build a hydrochemical hybrid code based on existing 
hydrodynamic codes by applying modules of chemical 
solvers and reaction networks at their appropriate time 
steps. Numerical method used in our hydrochemical hy¬ 
brid code is described in detail. The content is organized 
as follows. In Section[2l we describe basic equations that 
we shall solve, design of the code, and details of the nu¬ 
merical method. In Section [3l hydrodynamics module of 
the code is verified against hydrodynamic tests. In Sec¬ 
tion m chemical module of the cod e is verified against 
one dimensional PDR benchmarks inlRollig et alJ (1200'^ 
and multidimensional PDR model in iBisbas et al.l (1201211 . 
In Section [21 we present a hybrid case with both of the 
hydro-and chemical modules turned on, for a model of 
PDR on top of background thermal flows. In Section IH 
we discuss implications of the paper and future work. 


2 . NUMERICAL METHOD 
2.1. Basic Equations and Code Design 

The hydrodynamic equations for the conservation of 
mass, momentum, and total energy are written as 

^+V-(pn) = 0, (1) 


dv 


( 2 ) 


— + V ■ [(E + P)v] =T-A + pvg, (3) 

where p, v, P, g are the density, the velocity, the pres¬ 
sure, and the gravitational force, respectively. We in¬ 
clude the heating F and cooling A rates in the energy 
equation. We adopt the equation of state for ideal gas, 
so that the total energy density of the gas is expressed as 
E = P/{'y — 1) -I- p|t;p/2, where 7 is the ratio of specific 
heats. 

The chemical network includes the formation and de¬ 
struction reactions of all the considered chemical species. 
The abundance of each species i is controlled by a rate 
equation of the form 


duj 

dt 


I j k 


Ec.^+EE kijiUj 


(4) 


where Ui is the number density of species i, kj^i is the 
reaction rate coefficient for the reaction that form species 
i from the reaction of species j and fc, (^u is the local pho¬ 
todestruction rate coefficient for the ionization or disso¬ 
ciation of species i either by FUV photons or by cosmic 
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ray induced photons, producing species I as a result. The 
/iH denotes the mass per hydrogen nuclei whose typical 
value is ^ 1.15 mn in molecular cloud, where mn is the 
mass of atomic hydrogen, and the number density of hy¬ 
drogen nuclei tt-h is related with the gas density p by 
riK = p/ph- 

The operator-split method was adopted to solve hy- 
drodynamical and chemical evolutions self-consistently. 
The basic equations are solved by the hydrodynamics 
and chemistry modules. In our approach, evolution of 
the energy due to cooling and heating processes is sep¬ 
arately treated from evolution due to advection and ex¬ 
ternal forces. That is. Equation dJ]) is split into 

BE 

—-[{E + P)v]= pv ■ g, (5) 


and 




( 6 ) 


The hydrodynamics module solves Equations o, IH), 
and ([5]) . The number density of each species rii is treated 
as passively advected variables in this step. The chem¬ 
istry module solves radiative transfer to obtain local 
mean intensity of FUV radiation, and then updates the 
number density of chemical species and the internal en¬ 
ergy of the gas by solving Equations (ID) and ([6]). The 
algorithm adopted to couple the hydrodynamics and the 
chemistry modules is described in detail in Section [231 


is advantageous when flow in specific region of computa¬ 
tional domain needs to be solved with higher resolution 
than other regions. 

The hydrodynamics module includes gravitational 
forces due to self-gravity of the gas and a specified point 
mass. Gravitational potential for the self-gravity ^seif 
is obtained from Poisson’s equation 


VHself=^7TGp, ( 12 ) 

where G is the gravitational co nstant, and Equat ion (fT^ 
is solved by multigrid method ([Press et al.lll986H . To re¬ 
duce errors efficiently, iterative Poisson solver is run on 
the grid used for solving hydrodynamic equations and on 
coarser grids toget her. Boundary cond itions are deter¬ 
mined by following [Foster fc Bosa l|1996[l . The derivative 
d^/dr is set to zero on z-axis. On the outer boundary, 
^seif is calculated by Legendre polynomial expansion us¬ 
ing the density distribution in the computational domain. 
For the point mass its potential is given by 


^pm — 


GM, 


pm 


+ {z - Zpm) ’ 


(13) 


where Mpm and Zpm are the mass and z-coordinate of 
the point mass, respectively. The point mass is located 
on the z-axis by axisymmetry. Using these gravitational 
potentials, the gravitational force is obtained as 


g = -V($^e;/ +$pm)- (14) 


2.2. The Module of Hydrodynamics 

The module of hydrodynamics solves Equations o, 
and ([5|) . The change of total energy of the gas due to 
thermal process is solved in chemistry module. In cylin¬ 
drical coordinates system assuming axisymmetry, con¬ 
servative form of equations can be written as 


d rr 1 ^ 


where 


are conserved variables, 


and 


lively. 


dz^ 



(7) 

,Ef 

(8) 

,{E + P)Vrf 

(9) 

(E + P) 

(10) 

. z-direction, respec- 


5= 0, pgr 


P^l 


,P9z 


pV^Vr 


.pv-g 


( 11 ) 


are geometrical and physical source terms. These equa¬ 
tions are integrated by Godunov-ty pe finite volume 
scheme. An HLLC Riemann solver ([Toro et al.[ [199411 

is used to compute flux at the ce ll interfac e. _ Data 

reconstruction scheme proposed by [Mignonel (|2014[) is 
adopted to avoid large numerical errors near symmetry 
axis. Our hydrodynamics module has second-order ac¬ 
curacy in both space and time. Either uniform or non- 
uniform structured grid can be used. Non-uniform grid 


2.3. Chemistry Module 
2.3.1. Chemical Model and Thermal Evolution 

The module of chemistry solves a time-dependent set 
of equa tions for g a s-pha se che mistry based on t h e cod e 
used in [Lee et ah ([199611 and [Morata fc HerbstI (|2008il . 
The original FORTRAN source code was modified into 
a subroutine that could be repetitively called from the 
main routine of KM2. The evolution of chemistry was 
tested to be the same for different physical conditions of 
density and temperature in the integrated subroutine as 
in the original code. 

The chemical solver involves a system of highly non¬ 
linear algebraic equations of the type in Equation Q, 
one for each species. The publicly available Fortran sub- 
routine, DL SODEI^, which is based on the Gear method 
([Geail[197in . was adopted to find the time-dependent so¬ 
lutions. The Gear method is an implicit, linear multistep 
method that uses variable time steps and error control 
techniques to preserve the required accuracy during the 
integration. Equations for chemical evolution are solved 
on the same grid as for hydrodynamic equations. 

Arbitrary chemical reaction networks could be im¬ 
plemented through a Python routine in which differ¬ 
ent databases of chemical reactions could be read in. 
Recent databases for astrochemistry , such as UMISlQ 
( Woodall et al.[ [2007[: [McElrov et abi [2013[1 . and OSL0 
( Smith et al.[[2004 [Garrod et al.[[2008[l . include several 
thousand chemical reactions. The Python scripts read in 
reaction data set written in the data format of UMIST 

® http://computation.llnl.gov/casc/software.php 

^ http://www.udfa.net/ 

® http://www.physics.Ohio-state.edu/~{}eric/research.html 
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or OSU, and generate Fortran subroutines of computing 
time derivatives for all included chemical species. 

The change of the energy due to the heating and cool¬ 
ing in reactions is implicitly computed by backward Euler 
method. Update of the energy from time t to time t + At 
is expressed as 

E{t + At) = E(t) + At [r(t -h At) - A(t + At)]. (15) 


parallel cloud over all solid angle 

Xi(u2) = y^ Xo(q)exp(-7iyly(q)) —, (19) 

where q is the vector oriented to the cloud surface with 
respect to the position (r, z). Photoreaction rate at the 
position (r, z) is given by 


The cooling and heating rates are computed with chem¬ 
ical abundances at time t -I- At. The pressure and gas 
temperature are also updated using the updated energy. 
The gas temperature is given by 

MH[2n(H2)+n(H)+n(H+)] P 
n(H 2 ) -I- n(H) -|- n(H+) -|- n{e~) ksp" 

where ks is the Boltzmann’s constant. Hereafter, n(X) 
denotes the number density of species X. The mean 
molecular weight is obtained from number densities of 
H 2 , H, H+, and e“. 


2.3.2. Radiative Transfer and Photoreaction Rates 

In KM2, transfer of the FUV radiation is solved to de¬ 
termine rates of photoreactions such as photoionizations 
and photodissociations. The photoreaction rate due to 
FUV radiation is written as 

k= / a.^du, (17) 

Jut hv 

where is the cross section for reaction at the frequency 
i^Lyi is the Lyman limit frequency, vt is the threshold 
frequency for the reaction, h is the Planck’s constant, and 
Ju is the local mean intensity of FUV radiation. Since 
a frequency-dependent radiative transfer is computation¬ 
ally too expensive to be coupled with hydrochemical evo¬ 
lution at each time step, we adopt a frequency-integrated 
FUV intensi ty normalized with the Draine standard ra- 
diation field (lprainel|1978fl to obtain photoreaction rates. 
iRoberge et al.l ( 19911 1 showed that photoreaction rates of 
some species can be well fitted by function of visual ex¬ 
tinction by using frequency dependent radiative transfer. 
We assumed that this approximation hold for photoreac¬ 
tions in our chemical network. 

The rate of the ith photoreaction in a semi-infinite 
plane-parallel geometry is given by the form of 


kt{r,z) = aiXi{r,z). (20) 

Since photodissociation of CO and H 2 involves line ab¬ 
sorption, self- and mu tual shielding ef fects are important, 
and shielding factors (|Lee et al.lll996ll are included in the 
integration of Equation (IT^ for these photoreactions. 

2.3.3. Heating Processes 

FUV radiation is main heating source in PDRs. FUV 
radiation ejects electrons from dust grains due to pho¬ 
toelectric effect. The kinetic energy of ejected electrons 
goes into the thermal energy of gas through interactions 
with ambient hydrogen molecules. The photoelectric 
heating rate is (jBakes fc TielensI 1 199411 

Fpe = lO^^^eGo’^ ergs cm“^ s“^. (21) 

Here, Go is the mean intensity of FUV radi ation nor¬ 
malized with the Habing field (iHabind 1196811 . which is 
related with the intensity normalized with the Draine 
field by Gq = 1.7xo, and n = n(H) -|- 2n(H2) is the num¬ 
ber density of neutral hydrogen nuclei. The photoelectric 
heating efficiency e is written as 

4.87 X 10-2 

1 + 4 X 10-3(GoTi/2/„(e-))o.r3 
3.65 X 10-2(r/104)O-7 
1-^2 X 10-4(Gori/2/n(e-))’ ^ ’ 

Line absorption of a FUV photon will pump H 2 
molecules to a bound excited electronic state. About 
10 percent of the excited molecules drop back to the vi¬ 
brational continuum of the ground electronic state and 
dissociates, and other 90 percent drop back to an ex¬ 
cited vibrational state in the electronic ground state. 
Collisional de-excitation of FUV pumped H 2 can be im- 
portant heating source in d ense gas. Its heating rate is 
l|Hollenbach fc McKeelll97^ 


fc* = aiXoexp(-7iAv) (18) 

where ai is the photoreaction rate in the unshielded in¬ 
terstellar ultraviolet radiation field, xo is the scaling fac¬ 
tor of the radiation field at cloud surface with respect 
to Draine standard radiation field. Ay is the perpen¬ 
dicular visual extinction measured from cloud surface, 
and 7 i is the parameter representing attenuation proper¬ 
ties of the FUV. The visual extinction Ay is related to 
hydrogen column density Ah by Ay = 6.3 x 10~^^JVh 
(jWagenblast fc Hartauistl[l989H . 

For a cloud of arbitrary geometry, the photoreaction 
rates were computed in a similar manner as the three - 
dimensional PDR code developed bv iBisbas et al.l (1201211 . 
Local intensity of the FUV radiation within a frequency 
range contributing the ith photoreaction at position (r, z) 
is calculated by averaging that of a semi-infinite plane- 


Tj/y = 9Rrf(H2)—- , , eV cm'^ s-\ (23) 

1 -I- ncr(H2)/nH 

where i?ti(H 2 ) is photodissociation rate of molecular hy¬ 
drogen. The critical density is defined as 

ncr(H 2 ) = 10®r-i/2/{1.6xH exp[-(400/T)2] 

-fl.4a;2exp-[12000/(r-H 1200)]} cm-3(24) 

where xh = n(H)/nH and X 2 = n(H 2 )/nH are fractional 
abundance of atomic hydrogen and molecular hydrogen, 
respectively. When the density exceeds the critical den¬ 
sity, the excitation energy is effectively converted into 
heat via collisional de-excitation rather than into radia¬ 
tion. 

Photo destruction of atoms and/or molecules due to 
FUV radiation may also contribute to the heating of 
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gas. Photoionization of atomic carbon and photodisso¬ 
ciation of molecular hydrogen are considered in KM2. 
Kinetic energy of the ejected electron or H atom heats 
the gas. We adopted the value of 1.06 eV and 0.4 eV as 
released energies in a pho toionization of atomic carbon 
(|Meiierink fc: Spaaiisl 20051) and in a photodissociation of 
molecular hydrogen (iSpaanslllOOGl) . respectively. 

Cosmic rays become the main heating source where 
FUV radiation is strongly attenuated. Cosmic rays ionize 
molecules and atoms in the cloud. The kinetic energy of 
the ejected electrons is converted to thermal energy of 
the gas. Ionization of molecular hydrogen 

H 2 -b cr H -b H+ -f e“ (25) 

Hz -b cr H+ -b e“ (26) 

and ionization of atomic hydrogen 

H -b cr H+ -b e" (27) 

are considered as heating sources due to cosmic 
rays. With the cosmic ray ionization rates from 
iHollenbach fc McKed (|1989D , the cosmic ray heating rate 
is given by 

Ter = (0.952n(H2) -b 0.46n(H))()crA(5cr ergs cm“^ s“^, 

(28) 

where Ccr is the total rate for electron production from 
cosmic ray ionization and AQcr is the energy deposited 
as heat as a result of the ionizatio n. The values of Ar = 
5.0 X 10~^^ s“^ and AQcr = 20 eV (|Goldsmith fc Langed 
I1978D are adopted in the KM2. 

In interstellar clouds, molecular hydrogen is be¬ 
lieved to form mainly on dust grains. Formation rate 
of molecular hydrogen on dust grains is written as 
(|Tielens fc Hollenbachl[l98^ 


Table 1 

Atomic data for fine-structure transitions 


Species Transition A (/rm) A (s ’d 


C+. ^57 7 2.30 X 10-® 

C . 3Po - ®Pi 609.1 7.88 X 10"® 

®Pi - ®P2 370.4 2.65 X 10“’^ 

®Po - ®P2 230.3 1.81 X lO-i"' 

O . ®P2 - ®Pi 63.2 8.91 X 10“® 

3pi _ 3pp 145 5 1 75 X 10-5 

®P2 - ®P2 44.1 1.34 X 10"i° 


iTielens fc Hollenbachl [19851) . The cooling rate due to a 
transition from level m to level n is given by 


^xiymn) — Pm^mn^l^mn^esc(7'mn) 


Sxiymn) P{ymn) 
Sx {ymn) 


(32) 

where Um is the population density in level m, Amn is 
the spontaneous transition probability, ymn is the tran¬ 
sition frequency, (Sesci^mn) is the escape probability of 
photon at optical depth Tmm and P(vmn) is the Planck 
function at the background temperature of 2.7 K. The 
source function is written as 


Siymn^ 


( 9mnn 


(33) 


where c is the speed of light, gm and are the statistical 
weights of level m and level n, respectively. Level popula¬ 
tion is determined by solving the equations of statistical 
equilibrium 


Rran — E TljiRn 


n^m 


n^m 


(34) 


Rf(R 2 ) = 6 X 10-^^(T/300)°-^n(H)n5'(T), (29) 

and it depends on the sticking coefficient of atomic hy¬ 
drogen on dust grain 

S{T) = [1 -b 0.4(r -b Td)° ® -b 2 X 10“^T -b 8 x , 

(30) 

where Td is the d ust temperature d eterm ined by us¬ 
ing the method in IHollenbach et all (I1991D . This for¬ 
mation process yields 4.48 eV of binding energy, and 
4.2 eV and 0.2 eV of this energy are assumed to 
have been distributed into vibrational/rotational exci¬ 
tation of the hydrogen molecule and thermal energy of 
the hydrogen molecule leav ing the grain, respectively 
(IHollenbach fc McKeel 1197^ . Hence, heating rate due 
to formation of molecular hydrogen on grain surface is 

rgr = i?/(H2)|o.2+—- . . 1 eVcm-3 g-^. 

I 1 -bUcr(H2)/nH J 

(31) 


2.3.4. Cooling Processes 

Collisional excitation of atomic fine-structure transi¬ 
tion is an important source of radiative cooling in PDRs. 
The radiative cooling rates through the fine-structure 
lines of [O I], [Cl], and [C II] are calculated based on 
escape probability approximation (|de Jong et al.l 119801 : 


where 

Rmn — Amni^escij'mn^^^ ~b Qmn ) “1“ Cmn-i ^ 

Rmn — {dni 9 rn)-^nrnPesc{_'^nrn')Qnrn ~b Cmnt ijen n^5) 

Qmn — e ^i.^mn) 


Here, Cmn is the collision-induced transition prob¬ 
ability with atomic hydrogen and molecular hydro¬ 
gen. Atomic data for fine-structure transitions are 
taken fro m Leiden Atomic a nd Molecular Database 
fLAMDAPI lSchoier et al.l (|2005[ ) and summarized in Ta¬ 
ble [T] Let zj_ denote the depth normal to plane, the 
optical depth and the escape probability for a plane- 
parallel cloud with u niform density can be expressed as 
(|de Jong et al.l[l97^ 


i{z±) = 


An 


^T^^mn Jo 




nm{z'd)9n 


and 


Pescig'^ — 


1 — exp(— 3t) 


il Ek 

J Syd 

(36) 

(37) 


respectively, where Syd is the root mean square of the 
thermal and turbulent velocities. We use the escape 
probability averaging that of Equation (1571) over all solid 
angles. 


® http: //home. strw. leidenuniv. nl/ - {}moldata/ 
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Free electrons and protons are important coolants in a 
highly ionized region. When a free electron recombines 
with a proton, the emitted photon removes internal en- 
erg y of the gas. Cooling du e to recombination is given 
as (iHummer fc Seatonlll963ll 

A^ec = 03(1-09 + 0.158 X 10-'‘r)fcBTn(e“)n(H+) (38) 

where as = 2.7 x 10“^^ cm“^s“^ is the case B recom¬ 
bination coefficient for hydrogen. When a free electron 
is accelerated by a proton, its kinetic energy can be re¬ 
moved by emitting a pho ton. Cooling rat e due to free- 
free emission is given by (|Osterbrocklll989ll 

A// = 1.42 X 10-2^Ti/2g//n(e")n(H+) (39) 

where the Gaunt factor for free-free emission, g/j, is as¬ 
sumed to be 1.3. 

At the temperature of a few thousand Kelvin or higher, 
Lyo emission from atomic hydrogen and 6300 A line 
emission from atomic oxygen in the metastable level 
are effective cooling processes in atomi c gas. The co oling 
rate due to Lya radiation is given by (|SDitzei]ll978ll 

Alj^„ = 7.3 X 10-i®n(e")n(H)e-^i®’^°°/^. (40) 

The cooling rate due to [O I] 6300 A emission is given by 
((Tielens &: Hollenbachl[l985f) 

Aeaoo = 1-8 x 10-24n(O)[n(H) + n(H2)]e-22.800/T^ (4^^ 

Cool dust grains and carbon monoxide are important 
coolants in dense molecular regi on. The cooling rate 
of th e gas by cooler dust grain is (|Hollenbach fc McKed 

[T^ 


Arf = 1.2 X 


(—] 

\ioookJ 


1/2 


100 Ay 

^min / 


X [1 - 0.8exp(-75/T)] (T - Ta) ergs cm-^ s"\,42) 

where Oniin is the minimum size of grains. We adopt 
the value of 100 A for amin- Cooling due to CO ro¬ 
tational transitions is obt ained from tabulated cooling 
functions for T < 100 K bvINe ufeld et'ap l|1995ll and for 
T > 100 K bv iNeufeld fc KaufmanI ( 19931) . Optical depth 
averaged over all solid angles is used to obtain the cool¬ 
ing rate. For cooling due to vibrational CO transitions, 
collisional excitation rate t o v = 1 state is taken from 
iHollenbach fc McKed (|1989f) . The cooling rates from col¬ 
lisions with H and H 2 are 


Ago.vi = 3-0 X 10-^^A£;ior‘''^exp 


(^) 


-(^) 


3.43' 


X exp 

and 

A^i^j, = 4.3 X lO-^^AFlioTexp 


n(CO)n(H) 


3.14 X 10^ 


T 


X exp 


-3080\ 

) 


n(CO)n(H2), 


(43) 


0.333' 


(44) 


2.4. Time Integration 

In KM2, the modules of chemistry and hydrodynam¬ 
ics are coupled by operator splitting time integration 
method. In the hydrodynamics module, the time step 
of an update is determined by the Courant-Friedrichs- 
Lewy condition. 


‘^thyd = CcFL min 


Ar 


Cs -I- \Vr 


I\z 


(45) 


where Ccfl is a constant less than unity, Ar and Az are 
widths of cells along the r and z directions, and Cg is the 
sound speed. In the chemistry module, the number densi¬ 
ties of chemical species and thermal energy are updated 
in a more complicated manner. Equation Q is solved 
under the assumption of constant gas temperature and 
shielding factors for photoreactions. Spatial distributions 
of H 2 and CO affect their self-shielding. The timescale 
in which the self-shielding factors do not change signifi¬ 
cantly is given by 


^tsh (|rf„(H 2 )/dt|’|dn(CO)/dt|) ’ 

where Csh is a constant of order unity and Xq is the 
elemental abundance of carbon. The value of Cgh is typ¬ 
ically set to be from 0.1 to 1.0. Number densities of the 
chemical species are updated with the time step Xtsu, 
then thermal energy and temperature are updated using 
the cooling and heating rates computed with the updated 
chemical abundances. If the relative change of tempera¬ 
ture from the previous step is larger than the tolerance 
value Cth^ the chemical and thermal updates will be re¬ 
done with a smaller time step. The value of Cth is set 
to be less than 0.1 in the benchmarks. The thermal and 
chemical updates are sub-cycled relative to the single hy¬ 
drodynamic update, since timescales for changes of tem¬ 
perature and shielding factors are smaller than At^yd in 
general. 

Figure [1] shows the schematics of our time integration 
algorithm in block form. The procedure to update phys¬ 
ical quantities from time t to t + Athyd is summarized as 
follows. 


1. First, the time step Athyd for hydrodynamic up¬ 
date is computed from Equation (1451) . 

2. The physical quantities p, v, P, and E are updated 
with the time step Athyd by solving Equation ([J). 

3. The radiative transfer is solved to obtain FUV in¬ 
tensity and shielding factors for photoreactions. 

4. The time step for updating shielding factors, Atgh 
is computed from Equation (1461) . 

5. The time step for chemical and thermal update, 
Atch, is computed from Atch = At^h - Atch,tot, 
where At(.h,tot is the cumulative total of Atch- 

6. The number densities of included chemical species 
Ui are updated by solving Equation (|3]). The tem¬ 
perature and shielding factors are kept constant. 


respectively, where AEio = 3080 K fcs is the energy of 
transition from ?; = 1 to n = 0 states. 


7. The heating F and cooling A rates are computed 
with the updated chemical abundances. 
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Figure 1. Schematic diagram of the algorithm adopted. The 
hydrodynamics module and the chemistry module are boxed by the 
dashed lines, and the control loops are connected by solid arrows. 

8 . Then E and T are updated implicitly. If the rela¬ 
tive change of the temperature from previous step 
is larger than Cth^ go back to step Eland redo the 
following steps with smaller Atcft,- 

9. Steps E] through | 8 ] are repeated until Atch,tot be¬ 
comes equal to Atsh- 

10. Steps El through E] are then repeated until Atsh,tot 
becomes equal to Athyd, where Atsh,tot is the cu¬ 
mulative total of Atsh- 

Equations (j45ll and (j46ll are computed over all the com¬ 
putational cells, and their minimum values are adopted 
as time steps Athyd and Atgh- By contrast, the time 
step Atch is locally determined at each cell. It allows 
us to use large time step in cells where the gas is close 
to thermal equilibrium, and reduce computational time 
compared to using a global time step. 

3. TESTS FOR THE HYDRODYNAMICS 

The hydrodynamics module has been tested with some 
standard hydrodynamical problems with known analytic 
solutions such as the Sod shock tube and the Sedov solu¬ 
tion. The chemistry modu le was tur ned off during these 
tests. The Sod shock tube (lSodllI97§l is a Riemann prob¬ 
lem to test accuracy and ability of computational hy¬ 
drodynamics code in simulating compressible flow with 



Figure 2. The density (top), the pressure (middle), and the ve¬ 
locity (bottom) profiles for Sod shock tube test at t = 0.2. The 
diamonds represent numerical results. Analytic solutions are plot¬ 
ted by solid lines as references. 

shock wave. High pressure gas and low pressure gas were 
initially separated by contact discontinuity and at rest 
everywhere. The initial density and pressure were unity 
at z < 0.5, while 0.125 and 0.1 at z > 0.5. The ratio 
of specihc heats 7 was chosen to be 1.4. The size of the 
computational domain in z was unity and covered with 
150 cells. 

Figure E] shows profiles of the density, the pressure, 
and the velocity along z-axis at the time t = 0.2. Sod 
shock tube involves three type nonlinear waves, that is, 
shock wave, contact discontinuity, and rarefaction wave. 
The shock wave and the contact discontinuity propagate 
to the right, and the rarefaction wave propagates to the 
left. Positions of these three waves are identical to those 
of the analytic solution within resolution. The shock 
front is located at z = 0.85, and it is resolved with a 
few cells without post-shock numerical oscillation. The 
contact discontinuity can be seen at z = 0.68 only in 
the density prohle. The rarefaction front is located at 
z = 0.26. The density, the pressure, and the velocity 
behind it show good agreements with analytic solutions. 

As a multidimensional test, we show results of point¬ 
like explosion in homogeneous medium, which is known 
as the Sedov explosion problem. Sedov explosion involves 
self-similar propagation of a strong spherical shock wave 
through the background homog eneous mediu m. The ex¬ 
act analytic solution is given bv iSedovI (119590 . Our com¬ 
putational domain extended from 0 to 20 pc in z with 
1024 cells, and from 0 to 10 pc in r with 512 cells. The 
initial number density nn was set to be 1 cm“^. We de¬ 
posited a quantity of thermal energy E = 10®^ erg into a 
central small region whose radius is 0.4 pc to initiate the 
explosion. The ratio of specific heats 7 was chosen to be 
1.4. 

Figure El shows the number density field at t = 6 kyr. 
For presentation, the figure is symmetrized with respect 
to the z-axis. Shock wave spherically propagated, and 
good symmetry was maintained in the density field. We 
confirmed that the total energy and total mass of the sys¬ 
tem were conserved within errors of 2.5 x 10“^^ percent 
and 1.1 X 10“^^ percent through the simulation, respec¬ 
tively. Figure |4] shows density and velocity profiles along 
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Figure 3. The number density distribution at t = 6 kyr for Sedov 
explosion test. 



Figure 4. The number density and the velocity profiles for Sedov 
explosion test at t = 6 kyr as functions of the distance from the 
center of explosion. The dashed (green), the dotted (red), and 
the dash-dotted (purple) lines represent profiles along directions at 
angle 0, 45, and 90 degree with r-axis, respectively. The solid lines 
represent analytic solutions. 


different directions from the center of explosion with ex¬ 
act analytic solutions. Shock front is sharply resolved, 
and its position is identical to analytic solution. Al¬ 
though the velocity profile along z-direction slightly de¬ 
viates from the analytic solution, numerical results show 
good agreements with analytic solutions on the whole. 
The discrepancy between numerical and analytic solu¬ 
tions near z-axis is attributed to errors from geometrical 
source term ^ in Equation (0. Sedov explosion is the 
case that errors from geometrical source term ^ severely 
affect results, because high pressure gas was set at small 
central region as initial condition. In most other cases, 
errors from geometrical source term are negligible. 

4. TESTS EOR THE CHEMISTRY MODULE 

Similar benchmark tests have been conducted for 
the chemistry module with the hydrodynamics mod¬ 


Table 2 

Atomic initial abundances 


Element 

n(X)/nH 

H 

1 

He 

0.1 

C 

1 X 10-4 

0 

3 X 10-4 


ule turned of f. In this c o ntext, the PDR benchmarks 
presented in iRollig et al.l (|2007l) have been chosen to 
be the references. These tests mainly show how well 
the included chemical network based on the PDR mod¬ 
els and the specific cooling and heating processes are 
doing in comparison with other codes and their re¬ 
spective chemistry and physical processes designed to 
model PDRs. These tests (in the context of iRollig et al.l 
(|2007D 1 use a reduced chemical network only with some 
of the more fundamental molecules and reactions that 
had b een previous l y ado pted in the joint benchmark 
effort iR.ollig et al.l (j2007D by several different groups. 
In addition to these one dimensional PDR benchmark 
tests, applic a tion example of ‘^3D-PDR” presented by 
iBisbas et al.l (|2012D has been chosen for test of ability to 
solve multidimensional problems. Chemical and thermal 
structures of spherical cloud illuminated by EUV radia¬ 
tion has been solved in this test. 

4.1. The PDR benchmarks 

The reaction network a dopted exactly identical to the 
one in iR.ollig et al.l (120071) prepared for the benchmarking 
of PDR models includes the four most abundant elements 
(H, He, O, and C) for a total of 31 species. The reaction 
rates were taken from the UMIST99 database, with some 
corrections, for a total of 287 reactions. Table[2]gives the 
initial abundances of the species used, in which elemental 
cosmic abundances for H, He, C, an d O were ass u med. 

The eight test cases discussed in IR.ollig et al.l (j2007f ) 
were conducted, which covered number densities nn of 
10 ^ and 10® ® cm“®, and the impinging radiation fields 
with EUV intensities x of 10 10® in units of the 

Draine field. There were also two sets of models: the 
“E models”, in which a fixed gas temperature of 50 K 
and dust temperature of 20 K were adopted, and the “V 
models”, in which the gas and dust temperature along 
the PDR were self-consistently calculated. 

4.1.1. The F Models 

Figure [5] shows the comparison of results from the KM2 
code with ^^Lee96mod ’, a time-dependent code w hich 
uses the same treatment as iMorata fc HerbstI (j2008f ) and 
this paper, “UCL_PDR^’, another time-dependent code, 
“Leiden” and “Meudon ”, for the case of a fi xed gas tem¬ 
peratu re, r = 50 K, the ”F” Models, from IRollig et al.l 


The results are practically identical for most of the vi¬ 
sual extinction values in the four ”F” Models. When the 
cases of KM2 differ somewhat more from others, there is 
a wider spread of results among the models. At very low 
extinctions for the abundance of H 2 for the El model 
or at Av ^ 0.1 for the F4 model, the model predic¬ 
tions scatter while agree with one another reasonable 
well otherwise. For these cases, KM2 tends to fall within 
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Figure 5. Chem ical abundances as a function of visual extinction for the benchmark models with a fixed temperature of 50 K from 
IHoIlig et al.l II2007I1 and the KM2 code. The thick solid lines show the abundances of H (black), H2 (red), (green), C (orange), and CO 
(purple) from the KM2, and results from other models are labeled as ‘^LeeQdmod’ (dashed lines), ‘^Leiden!’ (dotted lines), ‘^UCL^PDR” 
(dot-dashed lines), “Meudon” (three-dots-dashed lines). Upper row) Model FI: density riH = 10® cm~®, FUV radiation field intensity 
X = 10 in units of the Draine field; Model F2: nn = 10® cm“®, x = 10®. Bottom row) Model F3: density njf = 10® ® cm“®, FUV radiation 
field intensity x = 10; Model F4; uh = 10® ® cm“®, x = 10®. 


the spread of abundances. The differences could be at¬ 
tributed to differences in the ways of calculating the H 2 
and CO photodissociation rates among the codes and not 
to the way that KM2 propagates the radiation from the 
outer to the inner layers. 

4.1.2. The V Models 

Figure [6] shows the comparison of the KM2 code for the 
case where the gas temperature is self-consistently cal¬ 
culated from the cooling and heatin g functions included 
in the codes. The cases taken from iRollig et al.l (|2007f) 
are “UCL_PDR”, ^Meiden” , ''''COSTAR” , and “Meudon”. 
The gas temperature will be different at each layer of the 
PDR. Figure [7] shows the distribution of the gas temper¬ 
ature as a function of the visual extinction into the cloud 
for the codes shown in Figure [51 

Figure [S] illustrates the models with variable temper¬ 
ature tend to show a significantly larger spread in the 
results among the codes. Still the KM2 in general agrees 
very well with the rest of the codes, for the ranges of Ay, 
and for the models VI to V3. They appear very simi¬ 
lar to those found in the cases in the ”F” models. On 
the other hand, the models vary greatly for some of the 
molecules (FI 2 , CO) in the V4 models, which have lower 
agreements among the participating codes in the range 
of Ay = 0.01 to 3. KM2 is not much worse than others 
in this regard. 

In the models VI to V3, the gas temperatures calcu¬ 
lated by the KM2 code are lower by a factor of a few for 
Ay < 1. This may be due to the larger cooling rate of 
atomic oxygen in KM2, whose collisi onal excitation rates 
of atomic oxygen were taken from lAbrahamsson et al.l 


(|2007D . These excitation rates are larger than values used 
in other co des previously by a factor of a few, which were 
taken from iLaunav fc Rouefl (1197711 . The cooling in O I 
is more effective than those of other codes. The V4 model 
shows again the largest difference in temperatures, and 
the temperature predicted by the KM2 code is within 
the spread of temperatures from other codes. These all 
indicate that performance based on the PDR specific 
pr oblem and the spe cific set of reaction rates adopted 
in IRollig et al.l (1200711 work well the overall design of the 
KM2. In next section, we discuss a sample model of its 
future use. 

4.2. The spherical cloud illuminated by a plane-parallel 
FUV radiation 

In this section, we present chemical and thermal struc¬ 
tures of a spherical cloud illuminated by a plane-parallel 
FUV radiation. The spherical cloud had uniform density 
of riH = 10^ cm“® and a radius of 5.15 pc. Our com¬ 
putational domain extended from 0 to 10.3 pc in z with 
2048 cells, and from 0 to 5.15 pc in r with 1024 cells. A 
center of the cloud was located at r = 0 pc and 2 ; = 5.15 
pc. The incident FUV radiation with intensity of x = 10 
entered from bottom boundary of the computational do¬ 
main. The density and intensity of the incident FUV 
radiation were same as model VI shown in Section [4.1.2l 
Chemical network and initial chemical abundances used 
in this test were also same as PDR benchmark tests. The 
evolution was solved for 30 Myr. It was long enough to 
achieve chemical and thermal steady state. 

Fig. |8]shows the temperature distribution of the cloud. 
For presentation, the figure is symmetrized with respect 
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Figure 6. The ”V” models of the chemical abundances as a fu nction of visual ext inction with a variable temperature solved through 
the heating and cooling processes for the benchmark models from fRollig et a n 1200711 and the KM2 code. The thick solid lines show the 
abundances of H (black), H 2 (red), C”*" (green), C (orange), and CO (purple) from the KM2, and results from other models are labeled 
as ‘^Lee96mod' (dashed lines), ‘^Leiden” (dotted lines), ^^UCL^PDR” (dot-dashed lines), “Meudon” (three-dots-dashed lines). Upper row) 
Model VI: density nn = 10^ cm“®, FUV radiation field intensity x = 10 in units of the Draine field; Model V2: n^i = 10^ cm“®, x = 10®. 
Bottom row) Model V3: density nn = 10® ® cm“®, FUV radiation field intensity x = 10; Model V4: nn = 10® ® cm“®, x = 10®. 



Figure 7. Temperature as a function of visual extinction for the 
benchmark models and the KM2 code (thick black line). Other 
codes are labeled as ^^COSTAR” (dashed lines), “Leiden^' (dot¬ 
ted lines), “UCMPDR” (dot-dashed lines), “Meudon” (three-dots- 
dashed lines). 

to the z-axis. The cloud surface of the bottom semi¬ 
sphere is heated to ~ 80 K by heating due to FUV ra¬ 
diation. As the FUV radiation is attenuated, its heating 
becomes inefficient. The temperature decreases to ~ 10 
K within thin region of 1 pc near the surface. The 


inner region of the upper semi-sphere has slightly higher 
temperature than outer region, because the efficiency of 
the cooling is lower there due to larger optical depth. The 
temperature profile along z-axis shows good agreement 
with results of model VI, the deviation between them is 
within a few K. Fig. [5] shows the chemical abundances 
as a function of visual extinction along z-axis and a di¬ 
agonal direction with 45 degrees respect to z-axis. The 
chemical abundances along z-axis show good agreements 
with those of model VI. On the other hand, intensity of 
FUV radiation decreases more steeply along the diago¬ 
nal direction. All transitions of H/H 2 , C^/C, and C/CO 
occur at smaller visual extinction than model VI. 

Our result s are consistent wit h those of “3D-PDR” 
presented in iBisbas et al.l (j2012[) . except the tempera¬ 
ture of surface region at the bottom semi-sphere is lower 
than ~ 170 K of ^‘3D-PDR”. ^^3D-PDR” is an extension 
of ^^UCL_PDR” for three dimensional models, therefore 
both codes produce identical results in one dimensional 
PDR benchmarks. As seen in Fig. [3 these codes tend to 
produce higher temperature than other PDR codes for 
models with FUV radiation intensity of y = 10. 

5. APPLICATIONS 

As a test example for the more general astrophysical 
problems, the hybrid hydrochemical evolution of a PDR 
is presented here for the intended capabilities of the KM2 
code. The initial density and intensity of FUV radia¬ 
tion field were set to be the same as the PDR bench¬ 
mark model V2, namely n = 10^ and x = 10^- The 
chemical network and initial chemical abundance were 
also the same as the benchmark models in Section |4l 
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Figure 8. The temperature distribution of spherical cloud illumi¬ 
nated by a plane-parallel FUV radiation. The dashed and dotted 
lines represent directions along which chemical abundances plotted 
in Fig. [9] 



of implemented chemical network affects hydrochemical 
evolution. 

Figure [To] shows the time evolution of the number den¬ 
sity, temperature, and velocity profiles. The incident 
FUV radiation enters from the left boundary and heats 
up the gas to ^ 4000 K at the surface. The expansion 
of the hot gas generates an evaporation flow that ac¬ 
celerates leftward at a velocity of ~ 6 km s“^. At the 
same time, a back-reaction of the photoevaporation flow 
drives a shock wave rightward and compresses the PDR. 
A dense compressed region with a temperature about 30 
K can be seen behind the shock front at t = 3 Myr. The 
shock fronts are located at z = 2.6 pc and z = 5.8 pc at 
t = 1 Myr and t = 3 Myr, respectively. From these val¬ 
ues, the propagation speed of the shock wave is 1.6 km 
s“^. The reference model with larger chemical network 
shows different evolution from PDR benchmark model 
V2. The photoevaporation flow in the reference model 
has higher temperature of ~ 1000 K and higher velocity 
of ~ 9 km s“^. The propagation speed of the shock wave 
is 2.0 km s“^. The stronger shock wave forms denser 
compressed region behind the shock front. 



0 1 2 3 , ,4 5 6 7 

z (pc) 


Figure 10. Number density, temperature, and velocity profiles as 
a function of z-coordinates in our hydrochemical simulations. The 
solid and the dashed lines represent results for PDR benchmark 
model V2 and for the reference model with larger chemical network, 
respectively. The dotted lines represent profiles at t = 0. 
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Figure 9. The chemical abundances as a function of visual ex¬ 
tinction along z-axis (dashed line in Fig. 11 and diagonal direction 
(dotted line in Fig. [Si . The dotted lines represent results of model 
VI in Section 14.1.21 


Our computational domain was one-dimensional and ex- 
tended from 0 to 7 pc in z with 1024 uniform cells. As 
a reference, simulation with larger chemical network has 
also been performed. For the reference model, 2799 re¬ 
actions composed of H, He, O, and C were taken from 
UMIST RATE12 database. The reaction rates for for¬ 
mation and photodissociation of H 2 were same as PDR 
benchmark tests. The chemical network for the reference 
model include 2801 reactions among 203 species in total. 
This reference model was intended to see how difference 


Figure [TT] shows the fractional abundances of H 2 , H, 
C’*', C, and CO at t = 1 Myr and t = 3 Myr. The molec¬ 
ular and atomic abundances clearly show the position of 
the shock front at both times, especially by the steep in¬ 
crease in abundance of CO. Since we assumed that the 
initial abundances were atomic, we can see the progres¬ 
sive increase of CO abundance in the inner shielded lay¬ 
ers of the cloud in time; and, similarly, the decrease of H 
due to conversion into H2. Figure [TT] also clearly shows 
how the increased density at the shock front enhances the 
formation of CO, which for t = 3 Myr is ~ 2 orders of 
magnitude higher than in the inner shielded layers. This 
result suggests that this contrast in abundance in CO 
may trace the position of the edge of the photoevapora¬ 
tion flow. The progression of the photo-evaporated flow 
is also traced, although not so dramatically, by the point 
of conversion of atomic into molecular gas, from H to 
H2. Figure [TT] also shows that the shock wave traversed 
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Figure 11. Chemical abundances at the time t = 1 Myr (top) 
and t = 3 Myr (bottom) as a function of z-coordinates in our 
hydrochemical simulation. The solid and the dashed lines represent 
results for PDR benchmark model V2 and for the reference model 
with larger chemical network, respectively. 

the entire computational domain before chemical equi¬ 
librium was reached in the uncompressed region. The 
time scale of hydrodynamic evolution is shorter than the 
time needed for chemical equilibrium. This affects the 
comparison with the models with variable temperature 
in Section Sm 

For comparison with results of model V2 shown in Fig¬ 
ure [HI we plot the fractional abundances of H 2 , H, C+, 
C, and CO at f = 1 Myr and t = 3 Myr as a function of 
visual extinction in Figure [121 The original V2 models 
were equivalent to the time when chemical equilibrium, 
in our case after 30 Myr. This difference in time accounts 
for the different abundances in the inner shielded layers. 
We can see in Figure 0 how the abundances at Ay 0.1 
are markedly smaller, by an order of magnitude or more, 
for H 2 , and C. After the passing of the photo-evaporating 
front, the density of this region is lower than that of the 
V2 case, which is reflected in lower shielding and a slower 
chemistry. On the other hand, at Ay > 3 the density is 
higher than model V2 by an order of magnitude due to 
shock compression and the effect of the higher gas den¬ 
sity on CO, C, and C+, both in terms of shielding and 
shorter chemical timescales can be seen clearly. The dis¬ 
sociation front of CO is located at Ay ~ 4, while it is 
located at Ay ~ 10 in model V2. The reference model 
shows different chemical evolution from PDR benchmark 
model V2 in some respects. One difference is faster con¬ 
version from C into CO in the shock compressed region 
because of higher density. Another difference is higher 
fractional abundance of H in the shock compressed re¬ 
gion. Chemical network for the reference model contains 
more reactions producing atomic hydrogen, such as colli- 
sional dissociation or cosmic ray dissociation of molecules 
including hydrogen. These reactions contributes higher 



Figure 12. Same as Figure [TT] except that the data are plotted 
as a function of visual extinction. 

fractional abundance of H. 

These hydrochemical hybrid simulations require about 
40 times longer computation time than pure chemical 
simulations without the evolution of the hydrodynam¬ 
ics. This longer computation time is mainly due to its 
shorter time steps determined by Equation (1451) . Com¬ 
putation times for chemical reactions, thermal processes, 
hydrodynamics, and radiative transfer account for 79.5, 
18.5, 1.2, and 0.1 percent of the total computation time, 
respectively. 

6 . DISCUSSION AND SUMMARY 

The basic concepts and algorithms that had been 
adopted into the hydrochemical hybrid code, KM2, have 
been presented in this paper. The code is modular, and 
consists of a hydrodynamics module and a chemistry 
module. These modules have been verified with hydro- 
dynamic test problems and PDR benchmark models that 
have been known in the respective fields for their perfor¬ 
mances. In particular, the method of time integration 
is advantageous in reducing the computation time that 
would otherwise be prohibitive for most problems. Be¬ 
cause the time steps for chemical and thermal updates 
are locally determined at each computational cell, and 
large time steps can be used in cells where the gas is 
close to thermal equilibrium. 

As an application example of the KM2 code, the first 
hydrochemical simulation of a photoevaporating PDR 
is presented in this paper. The photoevaporation flow 
changed the structure of PDR within a time scale shorter 
than the timescale needed to reach the chemical equilib¬ 
rium, that is the timescale that is shorter would happen 
first before the chemical equilibria are reached locally. 
Profiles of chemical abundances were different from mod¬ 
els without hydrodynamics. Moreover, photoevaporating 
PDRs with different chemical networks showed differ- 
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ent thermal and hydrodynamical evolutions. These sug¬ 
gest that hydrochemical hybrid models are important for 
studies of photoevaporation where multiple timescales 
are involved and physical processes are interacting with 
one another. In addition to the PDRs, protoplanetary 
disks irradiated by central star or FUV fields may be 
important targets affected by photoevaporation in which 
various timescales are competitive. Hydrochemical evo¬ 
lution models of photoevaporating PDRs and protoplan¬ 
etary disks will be investigated in separate papers. 

To apply the KM2 code to more astrophysical prob¬ 
lems, some extensions of its capabilities are planned as 
future work. In this paper, the chemical networks in¬ 
clude only gas-phase reactions, while grain-surface re¬ 
actions also play important role in chemical evolution in 
dense molecular clouds. Inclusion of the grain-surface re¬ 
actions will be desirable without the modification of the 
basic framework. The code can be parallelized by mes¬ 
sage passing interface (MPI) library for large scale sim¬ 
ulations. Reducing the computation time in the compu¬ 
tation of the chemical reactions is also essential for large 
scale simulations. As computations for chemical updates 
are done with local physical quantities, namely the den¬ 
sity, temperature, and FUV intensity, node communica¬ 
tion is not required in MPI parallelization of the chem¬ 
istry solver. Good scalability and improvement of per¬ 
formance can be expected once the MPI is implemented. 

Finally, this paper has demonstrated the essential 
implementation method for building a hybrid hydro¬ 
chemical code. This critical implementation allows the 
hydrodynamical and chemical evolution to be coupled 
by a straightforward operator splitting method in the 
KM2 code. The implementation method can be easily 
extended and applied to other hydrodynamical codes for 
their complete hybridization of the thermal evolution at 
each time step. 

The work has been funded and supported by the 
’’CHARMS” Initiative under the Theoretical Institute 
for Advanced Research in Astrophysics (TIARA) of 
Academia Sinica Institute of Astronomy and Astro¬ 
physics (ASIAA). Numerical simulations were carried 
out in part on PC clusters at the Center for Computa¬ 
tional Astrophysics, National Astronomical Observatory 
of Japan (NAOJ). 

REFERENCES 

Abel, N. R, Ferland, G. J., Shaw, G., &; van Hoof, P. A. M. 2005, 
ApJS, 161, 65 

Abrahamsson, E., Krems, R. V., &: Dalgarno, A. 2007, ApJ, 654, 
1171 

Aikawa, Y., Herbst, E., Roberts, H., & Caselli, P. 2005, ApJ, 620, 
330 

Asahina, Y., Ogawa, T., Kawashima, T., et al. 2014, ApJ, 789, 79 
Bakes, E. L. O., k Tielens, A. G. G. M. 1994, ApJ, 427, 822 
Bell, T. A., Viti, S., Williams, D. A., Crawford, I. A., k Price, 

R. J. 2005, MNRAS, 357, 961 

Bisbas, T. G., Bell, T. A., Viti, S., Yates, J., k Barlow, M. J. 
2012, MNRAS, 427, 2100 

Black, J. H., k van Dishoeck, E. F. 1987, ApJ, 322, 412 
Brinch, C., k Hogerheijde, M. R. 2010, A&A, 523, A25 


Cazaux, S., Tielens, A. G. G. M., Ceccarelli, C., et al. 2003, ApJ, 
593, L51 

de Jong, T., Boland, W., k Dalgarno, A. 1980, A&A, 91, 68 
de Jong, T., Dalgarno, A., k Chu, S.-I. 1975, ApJ, 199, 69 
Draine, B. T. 1978, ApJS, 36, 595 

Ferland, G. J., Korista, K. T., Verner, D. A., et al. 1998, PASP, 
110, 761 

Foster, P.N.,k Boss, A. P. 1996, ApJ, 468, 784 
Furuya, K., Aikawa, Y., Tomida, K., et al. 2012, ApJ, 758, 86 
Garrod, R. T., Weaver, S. L. W., k Herbst, E. 2008, ApJ, 682, 
283 

Gear, C. W. 1971, Numerical Initial Value Problems in Ordinary 
Differential Equations (Englewood Cliffs, NJ. Prentice Hall) 
Goldsmith, P. F., k hanger, W. D. 1978, ApJ, 222, 881 
Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 
439, 2386 

Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421 
Hogerheijde, M. R., k van der Tak, F. F. S. 2000, AkA^ 362, 697 
Hollenbach, D., k McKee, C. F. 1979, ApJS, 41, 555 
—. 1989, ApJ, 342, 306 

Hollenbach, D. J., Takahashi, T., k Tielens, A. G. G. M. 1991, 
ApJ, 377, 192 

Hummer, D. G., k Seaton, M. J. 1963, MNRAS, 125, 437 
Jansen, D. J., van Dishoeck, E. F., Black, J. H., Spaans, M., k 
Sosin, C. 1995, AkA, 302, 223 
Juvela, M. 1997, AkA, 322, 943 
Kamp, I., k Bertoldi, F. 2000, A&A, 353, 276 
Kamp, I., k van Zadelhoff, G.-J. 2001, A&A, 373, 641 
Launay, J. M., k Roueff, E. 1977, A&A, 56, 289 
Le Bourlot, J., Pineau Des Forets, G., Roueff, E., k Flower, D. R. 
1993, AkA, 267, 233 

Le Petit, F., Nehme, C., Le Bourlot, J., k Roueff, E. 2006, ApJS, 
164, 506 

Lee, H.-H., Herbst, E., Pineau des Forets, G., Roueff, E., k 
Le Bourlot, J. 1996, A&A, 311, 690 
McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&;A, 550, 
A36 

Meijerink, R., k Spaans, M. 2005, A&;A, 436, 397 
Mignone, A. 2014, Journal of Computational Physics, 270, 784 
Morata, O., k Herbst, E. 2008, MNRAS, 390, 1549 
Neufeld, D. A., k Kaufman, M. J. 1993, ApJ, 418, 263 
Neufeld, D. A., Lepp, S., k Melnick, G. J. 1995, ApJS, 100, 132 
Osterbrock, D. E. 1989, Astrophysics of gaseous nebulae and 
active galactic nuclei ( Mill Valley: University Science Books) 
Papadopoulos, P. P., Thi, W.-F., k Viti, S. 2002, ApJ, 579, 270 
Press, W. H., Flannery, B. P., k Teukolsky, S. A. 1986, 

Numerical recipes. The art of scientific computing 
Roberge, W. G., Jones, D., Lepp, S., k Dalgarno, A. 1991, ApJS, 
77, 287 

Rollig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187 
Sakai, N., Sakai, T., Hirota, T., k Yamamoto, S. 2008, ApJ, 672, 
371 

Schoier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., k 
Black, J. H. 2005, AkA, 432, 369 
Sedov, L. 1. 1959, Similarity and Dimensional Methods in 
Mechanics, ed. Sedov, L. 1. 

Shaw, G., Ferland, G. J., Abel, N. P., Stancil, P. C., k van Hoof, 
P. A. M. 2005, ApJ, 624, 794 

Smith, 1. W. M., Herbst, E., k Chang, Q. 2004, MNRAS, 350, 323 
Sod, G. A. 1978, Journal of Computational Physics, 27, 1 
Spaans, M. 1996, AkA, 307, 271 

Spitzer, L. 1978, Physical processes in the interstellar medium 
(New York Wiley-Interscience, 1978. 333 p.) 

Tafalla, M., k Santiago, J. 2004, A&A, 414, L53 
Tielens, A. G. G. M., k Hollenbach, D. 1985, ApJ, 291, 722 
Toro, E. F., Spruce, M., k Speares, W. 1994, Shock Waves, 4, 25 
van Dishoeck, E. F., k Black, J. H. 1988, ApJ, 334, 771 
van Weeren, R. J., Brinch, C., k Hogerheijde, M. R. 2009, A&A, 
497, 773 

Vanhala, H. A. T., k Cameron, A. G. W. 1998, ApJ, 508, 291 
Visser, R., Doty, S. D., k van Dishoeck, E. F. 2011, A&A, 534, 
A132 

Visser, R., van Dishoeck, E. F., Doty, S. D., k Dullemond, C. P. 
2009, AkA, 495, 881 

Wagenblast, R., k Hartquist, T. W. 1989, MNRAS, 237, 1019 
Wakelam, V., Herbst, E., Loison, J.-C., et al. 2012, ApJS, 199, 21 
Woodall, J., Agundez, M., Markwick-Kemper, A. J., k Millar, 

T. J. 2007, AkA, 466, 1197 


